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ABSTRACT 

Eikonal, Hamilton- Jacobi and Poisson equations can be used for economical nearest 
wall distance computation and modification. Economical computations may be 
especially useful for aeroelastic and adaptive grid problems for which the grid 
deforms, and the nearest wall distance needs to be repeatedly computed. 
Modifications are directed at remedying turbulence model defects. Eor complex grid 
structures, implementation of the Eikonal and Hamilton- Jacobi approaches is not 
straightforward. This prohibits their use in industrial CED solvers. However, both the 
Eikonal and Hamilton- Jacobi equations can be written in advection and advection- 
diffusion forms, respectively. These, like the Poisson’s Laplacian, are commonly 
occurring industrial CED solver elements. Use of the NASA CEL3D code to solve the 
Eikonal and Hamilton- Jacobi equations in advective-based forms is explored. The 
advection-based distance equations are found to have robust convergence. Geometries 
studied include single and two element airfoils, wing body and double delta 
configurations along with a complex electronics system. It is shown that for Eikonal 
accuracy, upwind metric differences are required. The Poisson approach is found 
effective and, since it does not require offset metric evaluations, easiest to implement. 
The sensitivity of flow solutions to wall distance assumptions is explored. Generally, 
results are not greatly affected by wall distance traits. 
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ABBREVIATIONS 


MF Marching front approach (used to speed up GGS) 

GGS Global Gauss-Seidel 

NNS Nearest normal surface search 

NSS Nearest surface search 

2D Two dimensional 

3D Three dimensional 
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1. INTRODUCTION 

Normal wall distances, d, are still a necessary parameter in a range of key turbulence 
models (see Fares and Schroder (2002)), Detached Eddy Simulation (DES) type 
approaches (Strelets (2001)) and also peripheral applications incorporating additional 
solution physics (Osher and Sethian (1988) and Tucker (2001)). Also, near wall iso- 
values can be used to convert, for example, tetrahedral meshes into hybrid ones. The 
lines of constant d form the basis of a body fitted grid (see Sethian (1999)). Far field d 
contours can also be used as a rapid means of evaluating computational interfaces on 
unstructured overset meshes having relative movements (Nakahasi and Togashi 
(2000)). Surprisingly, for highly optimised RANS/URANS (Unsteady Reynolds 
Averaged Navier-Stokes) solvers, the effort in calculating d can be a significant 
fraction of the total solution time. For example, even with Cray C90 class computers 
and time invariant meshes it can take 3 hours just to compute d (Wigton (1998)). 
Because of d evaluation expense, in some codes dangerous approximations are made. 
The following are given by Spalart (2000): 

a) computing distances down grid lines and not allowing for grid non-orthogonality; 

b) computing ‘J’ as the distance between a field point and the nearest wall grid point 
and 

c) in multiblock solutions ascertaining ‘J’ on a purely block wise basis ignoring the 
possibility that the nearest wall distance might be associated with another block. 

The latter can create large inaccuracies and also non-smooth d distributions that can 
inhibit convergence. In relation to point (c), Spalart (2000) notes for overset grids the 
situation can arise where the same point in space has different equations. 

Clearly the listed practices will give mostly inexact distances (or even multiple 

distances) d . However, the careful modification of d to some d can remedy 
turbulence model deficiencies or extend modeling potential (Secundov et al. (2001)). 

For example, d »d can remedy excessive sharp convex feature turbulence model 
influence (see Fares and Schroder (2002) and Tucker (2003)). For comers or 
bodies/surfaces in close proximity, there is an increased multiple surface turbulence 
damping effect (Mompean et al. (1996), Launder et al. (1975)). This effect can be 

modeled by setting d <d . 

Search procedures, integral approaches and differential equations can be used to find 
d. Crude searches require 0(n^nj operations where and n^, correspond to the 
number of surface and internal node points. Wigton (1998) and Boger (2001) present 
more efficient procedures needing 0(n^yln~) and 0(n^lognJ operations, 

respectively. Wigton’ s approach searches for nearest surfaces, i.e. it is a nearest 
surface search (NSS). Strictly, turbulence models are developed with the nearest 
normal surface (NNS) in mind. However, d fields for this are non-smooth. Integral d 
computation approaches are described in Boger (2001) and Launder et al. (1975). For 
complex geometries they are difficult to apply. Hence, the focus here is on differential 
equation based methods suitable for vector and parallel computers. 

There are several differential equation based methods. The level set method of Osher 
and Sethian (1988) solves a differential time (0 evolution equation 
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( 3<Z>/3? + |V(Z>| = rV^(Z> and F — > 0 ). The d field is evaluated in 0{n*) operations from 

the times at which the dependent variable reaches a set level. Motivated by solution 
expense, Sethian (1999) casts the level set approach into a boundary value framework. 
The Eikonal equation below is used. 

N=i (1.1) 

This can be solved in 0(n^ logn„) operations. The Equation (1.1) dependent variable, 

(j), models propagating front first arrival times. The right hand side implies the front 
has unit velocity, i.e. there is some velocity field with lUI = 1. Eares and Schroder 
(2002) essentially solve an Eikonal type equation for (j)~^ adding a modified Laplacian 

to give d . The Eikonal equation with a Laplacian, as below, is a Hamilton-Jacobi 
equation. 

|V<z)| = /(<zi)VV (1.2) 

Spalding (see Tucker (2003)) approximately reconstructs d from solution of a Poisson 
equation (V^^ = -l). The die/)) reconstruction involves an auxiliary analytically 
derived equation. An efficient Poisson equation solver typically scales as 
0{n^\ogn^). The Eikonal, unlike the Poisson equation (see Sethian (1999), Tsai 

(2002) and Tucker (2003)), is challenging to code. This is especially the case for 
unstructured grids (see Sethian and Vladimirsky (2002)). Hence its implementation in 
established industrial CFD solvers presents a significant time investment and hence 
cost. Therefore, here use of an Eikonal equation form that is amenable to general 
geometry CED code implementation is explored. The form is reminiscent of the 
Euler/Navier-Stokes equations (the key equations modeled in CFD solvers). Also, 
since the Poisson d method is easy to implement in such codes this approach is further 
studied. For the current work, as the base CFD solver for the wall distance studies, the 
NASA CFL3D program (see Krist et al. (1998)) is primarily used. 

Operation counts show differential approaches can be significantly faster than search 
procedures. However, for fixed mesh problems it is unusual for the search procedure 
to contribute significantly to the total CFD solution cost. Nonetheless, for the 
increasingly common flow solutions where the mesh deforms and sometimes also 
refines (perhaps with varying cell numbers) the repeated search cost becomes highly 
significant (see Boger (2001)). Then, differential approaches, which can make easy, 
safe use of previous d estimates, become especially attractive. The new d approaches 
presented here are best viewed in this moving mesh context. However, most of the 
examples chosen involve stationary meshes. These primarily test the robustness and 
convergence traits and so seem appropriate. The approaches considered are now 
outlined. 
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2. GOVERNING d EQUATIONS 

The Eikonal is an exact d equation. When defining the vector 


U = V^ 


( 2 . 1 ) 


it can be rewritten in the following advection analogous form 


\]»V<p = l. (2.2) 

The vector U corresponds to the front propagation velocity implied in the Eikonal 
equation and so again lUI = 1. As noted by Fares and Schroder (2002) and Tucker 
(2003), Laplacians can be used to modify d and so remedy turbulence model defects. 
Motivated by this, the following advection-diffusion d equation is also considered 

U*V</i = l + r(^)VV (2.3) 

Equation (2.3) when solved with (2.1) gives a Hamilton- Jacobi equation solution. 
The positive function F(^) is discussed later. Since (V^)^ = V*(^V^)-^V^^and 
U = , the following conservative form of (2.3) is also possible 

V*(<^) = l + r*VV (2.4) 

where 

r* = (r+<z)) (2.5) 

Along different lines, substitution of the Poisson based (J) distribution arising from 
solving Equation (2.3) with U = 0 and F = 1 into 


d=- 


T 

I 7 = 1 , 3 ' 


.^jj 


1 


I 


^jj 


+ 2 (/) 


( 2 . 6 ) 


also gives distances (see Tucker (2003)). The derivation of (2.6) assumes extensive 
(infinite) coordinates in the non-normal wall directions. Hence, unlike the Eikonal, d 
is only accurate close to walls. However, turbulence models only need d accurate 
close to walls. 


In summary. Equation (2.2) gives an Eikonal equation solution. However, to 
distinguish an Equation (2.2) solution from the direct solution of the Eikonal Equation 
(1.1), Equation (2.2) solutions are referred to as advective d equation solutions. 
Equation (2.3, 2.4) solutions are distinguished from direct Equation (1.2) Hamilton- 
Jacobi solutions by calling them advection-diffusion d equation solutions. Equation 
(2.6) solutions are referred to here as Poisson d equation solutions. 
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Laplacian form and role 

Near a fine convex feature (for example a wire), accurate distances are needed for 
theoretical correctness, so d-d . However, to prevent excessive far field influence 
d »d can be required (see Fares and Schroder (2002) and Tucker (2003)). Adjacent 
to a convex feature TV^d » 0 . Therefore, the positive Laplacian inclusion in (2.3) 
has the desired effect of enlarging/exaggerating d(=(p). Motivated by dimensional 
homogeneity, the need that as J — > 0 , J J but suggests 


T = £d (2.7) 

where £* is a constant. Clearly more ‘aggressive’ functions than (2.7) (e.g. 
r = £{-1 + 6"^) ) are possible but these are not explored. The Eikonal equation with 
is, as noted earlier, a Hamilton- Jacobi equation. 


At concave comers V^J<0, hence d«d. Therefore, with the Laplacian, the 
damping effects of ‘extra’ walls, discussed earlier, is naturally accounted for. At a J 
maximum, say for example at the center line of a channel, again <0. This is 
sensible (see Fares and Schroder (2002)), since again it will naturally model the 
damping effect of the two adjacent walls. Eikonal and hence Equation (2.3) distances 
will generally be discontinuous in gradient (e.g., this occurs where fronts effectively 
propagating from different walls collide). Therefore, the J term smoothing has the 
potential to enhance convergence. 
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3. NUMERICAL METHODS 

Here, the numerical methods to solve Equations (2. 2-2. 4) are described. The Poisson 
based method is not covered here. For numerical details on this method see Tucker 
(1998). Equations (2. 2-2. 4) can be solved on Cartesian, general curvilinear and 
unstructured grids. For curvilinear grids they must be transformed using the chain rule 
for differential calculus in, say, a ri) system. When solving in this, metric terms 
(M) such as ^^{=d^jdx) , and 77 ^ must be evaluated. This aspect will be 

discussed more later. Diffusive terms in Equations (2.3-2.4) are discretized using 
second order central differences. The advective derivatives in Equations (2.2-2.4) are 
discretized, just considering an v-coordinate direction for example, using the 
following order upwind type of approximation 


d(p 

dx 




(3.1) 


where 


Ki,i,A= 






Pm- P i 

^M 


(3.2) 


and 


= 1 for w > -u.^^ and u._^ > 0 else n._^ - 0 ; and 
=1 for -u.^^ >w._iand -u.^^ >0else = 0 . 


The propagation velocities in the above are simply evaluated using Equations (2.1), 
where, for example 


u 




Pi -P m 


Pm- P i 

^M 


(3.3) 


Setting r = 0 in (2.3) shows that the conservation form of the d equation has a 
Laplacian. Here the view is taken that the discretized conservation form of the 
equation is used (Equation (2.4)). However, for optimal solution of the advective d 
equation form (using an advancing front approach), the extra Laplacian must be 

ignored. The approximation of neglecting (^l^p could be considered as solving the 
non-conservation Equation (2.2) with the alternative approximation that 

M, (3.4) 

The above equation corresponds to using offset difference-based velocity 
components. These would correspond to those needed in the conservative Equation 
(2.4) form. Either viewpoint has no great accuracy implications. Where d needs to be 
accurate for turbulence modeling, Numerical tests confirm the validity of 

the above arguments. 
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Metric Discretization 

As shown in Tucker (2003), for strongly stretched grids the metric terms (M) must be 
carefully discretized. In CFL3D, for the i index ‘direction’ Mi_i and Mi+i metrics are 
first evaluated. The latter, for example, involves geometric data at the i and i+1 grid 
points. These Mi_i and Mi+i difference-based constructs are then, as with most CFD 
codes, averaged to give a single Mi value to be used in the discretised transport 
equations. However, as shown in Tucker (2003), for strongly expanding near wall 
grid spacings this results in d overestimations. The simple remedy used by Tucker 
(2003) is the following metric formulation 

M,. (3.5) 

The standard CFL3D implementation corresponds to = 1/2 in the above. 

Stabilizing Measures 

To ensure stable solution of the transport based d equations, velocity clipping and 
diagonal dominance enhancement have been used. Both use the observation that, in 
2D, the exact U field should satisfy 

K=\ul+vl-\\ = 0 (3.6) 

where u and v are velocity components in the x and y directions, respectively. 
Therefore, to improve diagonal dominance, and Ru^i are added to both the 
discretized equation matrix diagonals and source term, respectively. Based on 
Equation (3.6), the following velocity clipping is also used 

|m,,|<1, |v,.|<1 (3.7) 

The Eikonal equation does not permit a backwards front movement (see Sethain 
(1999)). The implication of this restriction is that if in (2.3), for example, TV^0<-\ 
a theoretical violation occurs (the sign of the Equation (2.3) right hand side gives the 
front propagation direction). Hence, when solving the advection-diffusion form of the 
d equation, the Laplacian of Equation (2.3) (L) is modified to 

L = max[-C,rVVj (3.8) 

where 0<C<1. Here, C = 1 is used. This is the limit for theoretical correctness. 
Since anti-diffusion is associated with instability, use of (3.8) should be viewed as a 
stability measure. 

As a further stability measure, under-relaxation is used either through the following 
^ = {I - -h (where a is an under-relaxation factor and the ‘new’ and ‘old’ 
superscripts indicate iterative states), or through the use of a pseudo time term. 

When solving the Eikonal equation, evaluation of the discretized matrix coefficients 
for each grid point needs around 160 operations. The advective Eikonal equation form 
(Equation (3.6)) needs 90 operations, but then iteration is also required (i.e. Equations 
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(2.1) and (2.3) must be contained in some iterative loop). Therefore, the key 
efficiency issue is the number of iterations N typically required. This is discussed 
later. 

Simultaneous equation solution 

Here, either essentially crude global Gauss-Seidel (GGS) type iterations or a marching 
front (MF) approach is used. The essence of the MF approach is illustrated using Fig. 
1. In this, a ‘wire’ is represented on a Cartesian grid using a single node point. Points 
surrounding this are considered to be either seed, trial, or those to be later solved. 
These points are labelled through the marker variable ns with the following respective 
values 2, 1 and 0. To start, a seed (or multiple seeds) point is defined. The advective d 
equation is applied N times (until convergence) to assign this point a d value. As 
noted earlier, the N iterations involve repeated application of the discretized d 
equation followed by substitution of the result in Equation (2.1) to find U. The seed is 
represented in Frame (1) by the closed symbol. The advective d equation is then used 
with Equation (2.1) for all immediate seed neighbors (trial points). These ns = \ 
points are shown as open symbols. As shown in Frame (II), the trial point with the 
minimum d is taken as the next seed. Now, for this point, ns = 2. Then, as shown in 
Frame (III), d for all immediate neighbors to the new seed (for which ns 2) are 
calculated. Hence the next seed point location is found. Frames (IV) and (V) show 
subsequent development stages for a clockwise moving circular front. In summary the 
procedure is as follows: 

a) Define a seed point; 

b) Iteratively solve for d at immediate neighbors to the seed (trial points) and 

c) Make the trial point with the minimum d the next seed and return to (a). 

This procedure is continued until ns -2 for all points in Q or a smaller domain/extent 

Q identified by the Fig. 1 dotted line (many turbulence models and DES only need d 
to a maximum of about 1/3'^^* the boundary layer thickness). In a practical system, all 
surface adjacent points can be taken as seeds, and what is called an ‘active front’ is 
produced. For best efficiency a heap-sort procedure is required. This procedure has 
not been incorporated as part of this work. When a Laplacian is included, the MF 
approach can only be used to obtain an initial guess. Cleary, based on operation 
counts, A < 2 is needed for the advective d equation to match the Eikonal’s efficiency. 

Initial/Starting Conditions 

For the advection analogous d equations, i < ^ < oo in is found to be adequate. 
However, for MF solutions, (/> ^ o<> is used. When using the MF approach to give a 
DES distance field, (j) = Cdes^i is used where A, = max(Ax,Ay,Az) is the DES ‘filter’ 
control. The d computation naturally stops when d - Cdes^u For the Poisson 
approach, (p-Ois adequate. 

Boundary Conditions 

Boundary conditions on the domain boundaries are now described. At solid walls the 
following Dirichlet condition is applied 


^ = 0 


(3.9) 
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At flow/far field boundaries 



dn 


(3.10) 


can be used, where n is the boundary normal co-ordinate. However, if Q is 
sufficiently large, (3.9) makes a stable, most computationally economical far field 
boundary condition. It is especially preferred for the Poisson method, for which it 
gives much faster convergence. 



4. DISCUSSION OF RESULTS 

The following geometries are considered: 

(a) Flat plate; 

(b) Electronics system; 

(c) Single element airfoils (NACA 4412 and 0012); 

(d) Wing body; 

(e) Wing flap and 

(f) Double delta. 

For cases (c-f) where d ‘error’ values are given, these are relative to NSS distances. 
Of course the NSS cannot be regarded as exact wall distances but they do serve as a 
useful reference point. Negative errors correspond to a over-prediction. Unless 
otherwise stated, flow solutions use the Spalart and Allmaras (1994) (SA) turbulence 
model. This probably has the most extreme d requirements of most standard 
turbulence models. It needs d accurate for around l/3‘^‘' the boundary layer thickness. 
Hence it is a good candidate for d modeling tests. Quoted CPU times are for a 300 
MHz SGI Origin. 


Flat Plate( Case (a)) 

To test advective d equation solution method performance, a flat plate is considered. 
Other d methods are not considered for this case. The plate is positioned at y = 0 in a 
2D square domain Q. with sides of unit length. Equation (3.10) is used at the ‘far- 
field’ boundaries. Essentially, the two Fig. 2 grids are considered. The mildly 
sinusoidal function, distorted Frame (a) mesh is used as a basis for exploring the 
effect of grid expansion on accuracy. To do this the y direction grid positions are 
strongly expanded by raising them to some power k. The revised grid locations are at 
y^, where y is for the initial grid. In Frame (a) k = 1 . The Frame (b) grid also uses k = 
1. However, the frequency and amplitude of the sinusoidal disturbance is much 
higher. This grid is intended to severely test the accuracy and robustness of the 
advective d equation solution process. Table 1 gives % d errors for different k with 
offset (Equation (3.5)) and centered metrics. The parameter in the second column (Oy) 
is a geometric grid expansion factor estimate. 


k 

Oy 

% d Error 

Offset 

metrics 

Centered 

metrics 

1.1 

1.1 

0.48 

1.53 

1.5 

1.7 

0.99 

8.25 

2 

2.6 

3.15 

15.0 


Table 1. % d errors with offset 
and centered metrics. 


Clearly, with centered differences and larger k values, serious errors arise. This is 
perhaps not surprising. As noted by Roache (1976), for accurate flow solutions Oy < 
1.3 is needed. Also, the front propagation nature of the advective d equations makes 
errors additive. Eurthermore, centered differences are inconsistent with a propagating 
front problem. Eig. 3 plots predicted against actual distances (y) at x=0.5. Erame (a) 
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is for the Fig. 2a grid with k=2.6. Frame (b) is for the Fig. 2b grid. The full lines give 
the exact solutions. The circles are for offset metrics and the triangles for centered 
metrics. The Frame (a) d over-prediction arising from centered metrics is clear and 
consistent with the Eikonal equation observations of Tucker (2003). Frame (b) shows 
that (with offset metrics) even with an extremely distorted grid, reasonable distances 
can be obtained. 

Fig. 4 shows strongly distorted grid (Fig. 2b) velocity data. Frames (a) and (b) give 
lUI contours and U vectors, respectively. Considering the poor grid form, the velocity 
magnitudes and direction are surprisingly good. Fig. 5 gives a plot of log (Ru) (see 
Equation (3.6)) against GGS iterations for the strongly distorted grid. An under- 
relaxation factor, a= 0.7, is used. Through its monotonic nature, the plot reflects the 
stable convergence. 

Complex electronics system (Case (b)) 

Next, the Eig. 6 geometry is considered. Eor this case Poisson, NNS, Eikonal, 
advective and advection-diffusion d solutions are made. Eig. 7 gives v-y and y-z plane 
views of the Cartesian blocked off grid used as part of Case (b) studies. Initially some 
convergence studies are made using a coarsened Eig. 7 grid structure with just 20,000 
points. Eig. 8 uses a histogram to indicate advective d equation ME convergence. It 
shows the percentage of grid points needing various N values for an under-relaxation 
parameter (a) of 0.6. Generally A = 4. Therefore, since the advective d equation 
operation count is 45% less than the Eikonal, the former has roughly twice the 
computational cost. However, clearly, for non- stationary grids N will greatly drop. 
Also, as previously noted, for unstructured grids and industrial CED solvers the 
advective d equation is a more practical option than the Eikonal. Eig. 9 plots Ru 
against N for GGS with a =0.2, 0.6 and 0.8. Clearly, for GGS, a =0.8 gives best 
convergence. This high value reflects the advective d procedure’s convergence 
robustness. Table 2 gives ME to GGS total iteration ratios and ME scheme N values 
for different a. ME convergence is taken to occur when Ru < 0.001. Eor GGS this 
convergence level could not be achieved, Ru being around 30% higher. 


a 

MF/GGS 

7V(MF) 

0.2 

0.14 

9 

0.4 

0.17 

5 

0.6 

0.26 

4 

0.8 

0.47 

5 


Table 2. Iteration ratios and N. 


Table 2 shows the ME approach significantly reduces the number of iterations at each 
grid point. Also, pleasingly, N is not an excessively strong function of a. 

The Laplacian in the advection-diffusion form of the d equation (Equation (2.3)) had 
little convergence influence. Even with 10, the 6tr= 0.6 convergence is virtually 
identical. It can therefore perhaps be concluded that the advective d equation is quite 
stable, not needing any additional stabilizing diffusion. 

Eig. 10 gives circa mid v-y plane U and d error data for an advective d equation 
solution. Correctly, the Erame (a) vectors are surface normal. In corners, shock 
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analogous features can be found. At ‘convex’ corners these are the expansion type. 
Correctly, the frame (b) lUI contours are close to unity having an average value of 
1.004 (using a 0.001 % error in Ru as a basis for convergence). Frame (c) gives the % 
difference in d relative to a NNS field. The major discrepancies are seen to be around 
the ‘shock’ features and right hand side flow outlets. 

Fig. 11 gives x-y plane d contours. Frames (a-e) are for the NNS, Poisson, Eikonal, 
advective, and advection-diffusion solutions, respectively. The NNS contours are, as 
they should be, discontinuous. When using the CFL3D NSS procedure, sharp convex 
feature distances become similar to those for the Fig. 11 Eikonal/advective d equation 
results. Eor illustrative purposes, the advection-diffusion d solution (Erame (e)) uses 
an extremely large 6*= 10 value in Equation (2.7). Such a large value would (unless a 
function with a more rapid near wall decay than (2.7) is used) contaminate d accuracy 
and hence could not be used as part of a practical CED solution. The Poisson result 
naturally gives strong comer contour rounding. Correctly, the Eikonal and advective d 
equation solutions are virtually identical. Their distance traits appear as a compromise 
between those of the NNS and Poisson methods. Over the complete domain the 
average d deviation (relative to the search) for the Eikonal, Poisson, and advective 
equation approaches is 4.4, 17.1, and 4.8%, respectively. These differences appear 
significant. However, in the regions of greatest importance, near walls and away from 
corners, all the approaches give very similar distances. To see this, it is worth 
comparing Eikonal d solutions with the zonal k-s/k-l Poisson based d predictions of 
Tucker and Pan (2002). In this work, velocities and turbulence intensities are 
compared with measurements over six profiles. The average solution difference over 
these profiles, for the two d approaches, is just 2.8%. The discrepancy between the 
flow solutions and measurements is around 25%. Pig. 12 gives zoomed in views of d 
contours around regions A and B of the Erame (e) £= 10 solution. The Pig. 12 views 
show that the Laplacian has the effect of exaggerating d around sharp convex features 
(contour lines bias themselves to corners). Around concave features the reverse 
occurs. As noted earlier, both traits can be beneficial, remedying turbulence model 
defects. 

Pig. 13 shows the RANS/LES interface for a DBS solution. This was computed using 
the advective d equation and ME approach. To bring the interface away from the 
surface, so it can be seen more clearly, the Pig. 7 cell density is halved and Cdes = 1-3. 
This is twice the standard Cdes value. The advective approach with the ME solution 
naturally gives a DBS ‘ J’ field. It avoids wasted computational effort by terminating 
at the RANS/LES interface. The non-uniform nature of the interface reflects a DBS 
complex grid deficiency. DBS implies the interface location is at J = Cdes 
max(Ax,Ay,Az). Prom studying the Pig. 7 grid, it is clear this interface will be non- 
smooth. 

Pig. 14 plots d against y for the Poisson, search, and advective d solutions. The lines 
are located around the mid v and z locations of the system. The full, long and short 
dashed lines correspond to the search, advective, and Poisson approaches, 
respectively. As can be seen, away from walls the Poisson distances differ 
significantly from those for the other approaches. 
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Single element airfoils (Cases (c)) 

In this section, results for both NACA 4412 and 0012 airfoils are considered. The 
Poisson, advective J, NSS, and NNS procedures are compared. The NACA 4412 
angle of attack is fixed at 13.87^, M = 0.2 and. Re = 1.52 x 10^ based on chord. For the 
NACA 0012, the mean angle of attack is 4.86^, M = 0.6 and Re = 4.8 x 10^. The wing 
has a forced pitching oscillation of 50.32 Hz with an amplitude of 2.44^ (AGARD 
Case 3). Wall distances are computed using modified CFL3D solver elements. These 
modifications include use of upwind metric differences. For the NACA 4412 and 
0012, comparisons are made with the measurements of Coles and Wadock (1979), 
and Landon (1982), respectively. 

Fig. 15 shows the single block ‘C’ and 3 zone overset grids used for the NACA 4412 
studies. Both involve around 25,000 grid nodes. Fig. 16 shows circa 33,000 cell 
NACA 0012 ‘C’ grids at around the minimum (2.6^) and maximum (7.2°) pitch 
angles. With the grid re-distribution used (see Bartels (2000)) the near wall grid tends 
to move with the surface. Hence, for this simple case, the near surface d values should 
alter little. 

First, we consider the NACA 4412 case. Fig. 17 shows d contours representing in the 
following respective frames the: (a) NSS; (b) NNS; (c) advective equation; (d) 
Poisson and (e) overset grid Poisson results. Frames (a)-(d) represent results on the 
‘C’ grid, and Frame (e) represents results on the 3-zone overset grid. Overset Poisson 
solutions are made for the complete domain Q , and also just for the surface mapping 

block region (Q). Obviously, the latter solution, shown in Frame (e), is most 

computationally efficient. At the Q boundaries the Dirichlet condition ^ = 0 is used. 
Perhaps the most strikingly different d contours are those representing the nearest 
NNS shown in Frame (b). 

Fig. 18 gives y'^<400, d error histograms. The error is defined using the equation 
below 


Error = ) 

^NNS 


(4.1) 


where dNNS is the distance from the nearest surface search. The d without the subscript 
in Equation (4.1) is the distance field for which the error is being evaluated. Frame (a) 
is for the advective d equation with upwind metrics. The average ‘error’ is 0.5%. 
Without upwind metrics d is overestimated by an average of 6%. Frame (b) is for the 
overset grid Poisson solution using Q . Although not shown by the error histograms, 
the average d error for the Poisson is significantly higher than that for the advective 
equation. The ‘C’ grid Poisson error histogram looks very similar to that for Frame 
(b). 

Fig. 19 gives a zoomed in view of trailing edge region Poisson d ‘error’ contours. As 
can be seen, a key d overestimation zone is the sharp convex trailing edge geometry 
region. This overestimation zone is potentially desirable. As can be seen later, 
distance overestimation can remedy a turbulence model defect. 
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Fig. 20 gives U vectors for the advective d equation. Frames (a) and (b) are U vectors 
for the ‘C’ and overset grids, respectively. The vectors correctly show surface normal 
U components. The magnitude of these is close to unity. 

Table 3 gives NACA 4412 lift coefficients (Cl) and drag coefficients (Cl>) for the 
advective, NSS, NNS, and Poisson methods. 



Cl 

Cd 

Advective 

1.704 

0.03466 

NSS 

1.698 

0.03496 

NNS 

1.712 

0.03514 

Poisson 

1.713 

0.03543 

Table 3. NACA 441 

.2 Cl and 


Cl)’s for different d fields. 

Some minor differences are evident. Their nature is partly clarified in Fig. 21, giving 
surface Cp (pressure coefficient) curves. The full, dashed and dot-dashed lines are for 
the Poisson, NSS, and NNS approaches, respectively. The symbols give measured 
data. Differences can only clearly be seen when zooming in close to the trailing edge 
region. The advective d equation results are not shown. However, they do not fall 
outside the curve envelope for the other methods shown in Fig. 21. 

Fig. 22 plots top airfoil surface streamwise mean fluid velocity against the cross- 
stream coordinate. The plot is made near the trailing edge at x/c = 0.953 (c is the 
airfoil chord and x/c = 0 corresponds to the leading edge). This is the only region 
where any velocity discrepancies can be discerned. The velocity is normalized by the 
free stream value Uq. The line styles are the same as in Fig. 21. Velocity differences 
for the different d fields are small, confined to the trailing edge region and of no great 
consequence. 

The SA model’s turbulence destruction terms can be excessively active around sharp 
convex features (see Tucker (2003)). Fig. 23 gives SA trailing edge region turbulent 
viscosity (//t) contours for the NSS (Frame (a)), NNS (Frame (b)), and Poisson (Frame 
(c)) approaches. The large NNS distances in the region x/c > 1 greatly reduce the 
modeled turbulence destruction. As noted earlier the Poisson method also has this 
trait. The reduced turbulence destruction for the NNS and Poisson approaches relative 
to the NSS is evident in Fig. 23. Fig. 24 compares Poisson predicted trailing edge 
Reynolds stress contours with those inferred from measurements. Frames (a), (b), and 

(c) give the following respective components: uujul, uvjul and Wjul. For 

predictions these values are approximated from where S[j is strain 

rate and p fluid density. Contour levels have been chosen to match those identified for 
the hot wire measurements. As might be expected, due to the simpler nature of the SA 
model (for which the predicted stresses have very limited meaning) and vagaries of 
contour plotting the correspondence between the predictions and measurements is not 
high. However, both sets of results suggest higher turbulence energy zones for 

and 
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Fig. 25 plots uujul against y/c at x4: = 1.1. The full, dot-dashed, and dashed lines 
are for the Poisson, NNS, and NSS approaches, respectively. Symbols give measured 
values. Clearly, the NNS and Poisson approaches, having larger x/c> I d fields and 

hence less modeled destruction, give larger uujul values. The advective d results, 
not shown in the graph, are similar to those for the NSS. Table 4 gives 
integrated/averaged % normal stress intensity values, for -0.05 < y/c < 

0.028. As suggested by Fig. 25, the Poisson and NNS approaches yield a higher peak 
that is somewhat closer to the measurements. Again, as noted for Fig. 24, SA has 

limited accuracy when predicting uujul . Hence, the above results should be partly 
viewed as giving quantative information on the sensitivity of turbulence statistics to 
wall distance assumptions. 



(“/ ) % 

Measured 

3.10 

Poisson 

2.90 

NNS 

2.89 

Advective 

2.77 

NSS 

2.77 


Table 4. Average x/c =1.1 
normal intensity values. 


Clearly, for the current case, turbulence intensities downstream of the airfoil have 
insignificant influence on the parameters of interest in an aeronautical design context. 
However, for multiple element airfoils it is not inconceivable that more significant 
solution differences could arise. These can be caused by a peak in modeled turbulence 
energy convecting close to the center of the leading edge of a downstream element. 
Small changes in the turbulence energy peak’s position could give rise to different 
solutions (Spalart (2002)). 

Next, we present results for the pitching NACA 0012 case. Fig. 26 compares NACA 
0012 lift and moment coefficient (Cm) measurements, represented by symbols, with 
predictions. The full lines give an NSS d solution with a non-deforming grid that 
rotates with the airfoil. The long dashed and dotted lines are for deforming grid NSS 
and Poisson d distribution solutions, respectively. Importantly, the NSS solution uses 
a time invariant d distribution taken from the non-deforming grid case. Frames (a) and 
(b) give lift and moment coefficient data, respectively. For time accurate solutions, at 
least 10 sub-iterations per time step are required. Streamwise Courant numbers are 
around unity. The similarity of the lines suggests, as might be expected for this simple 
geometry, that the scheme used for the grid redistribution keeps the surface grid form 
quasi-constant with time. Hence it is very similar to that for the fixed mesh case. 
Therefore, the major differences between the fixed and deforming mesh results can 
most probably be attributed to discretization errors arising through the mesh 
movement and not wall distance differences. As might be expected, the Poisson 
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equation approach needs little convergence effort. However, this is to a large extent 
due to the relatively low levels of grid deformation for this case. Nonetheless, the 
Poisson approach is found to require a factor of around 2.25 less computing time than 
the search. Also, GGS solver iterations comprise about a factor of 1/100 of the total 
Poisson cost. Therefore, with large grid movements the Poisson method is likely to be 
very efficient. 

Wing body (Case (d)) 

For this wing body case the angle of attack is 2.87^, M = 0.802 and Re = 13.1 x 10^ 
based on the wing chord. The single block grid shown in Fig. 27 involves around 0.9 
million cells. For this case the Poisson, advective d, and NSS d procedures are used. 
Fig. 28 gives the < 400-region ‘error’ histogram for the advective d and Poisson 
methods. The average ‘errors’ are 3.13% and 3.14%, respectively. The Poisson’s 
tendency to over predict d is evident in the histograms. For this more complex 
geometry case, unlike the NAG A 4412 results, the average differences between lift 
and drag coefficients for the different d fields is lower. This can be seen in Table 5, 
giving Cl and Cd values for the different d fields. 



Cl 

Cd 

Advective 

0.6633 

0.04839 

NSS 

0.6632 

0.04855 

Poisson 

0.6635 

0.04842 


Table 5. Case (d) Cl and 
Cl)’s for different d fields. 


Wing Flap (Case (e)) 

The Case (e) wing and flap angles are 2^ and 40^, respectively. The wing-flap gap is 
0.6% of the wing chord (c). The Reynolds number, based on c, is 23 x 10^ and 
M=0.18. Part of the 10-block, approximately 0.9 million-cell grid is shown in Fig. 29. 
For this case the Poisson, NSS, and advective d procedures are tested. 

In Fig. 30, d contours for the NSS (Frame (a)), advective (Frame (b)), and Poisson 
(Frame (c)) methods are qualitatively compared. As can be seen, near walls the 
contours are similar. For the Poisson method, away from walls there is some unsightly 
d behaviour. This, as noted in Tucker (2003), is due to significant levels of block 
interface grid non-orthogonality (the CFL3D solver assumes grid orthogonality) and 
the use of backward differences in Equation (2.6). The Poisson equation’s trait of 
overestimating trailing edge d values is evident. Fig. 31 plots the d equation residual 
log against the number of GGS type iterations. As can be seen, the residual can be 
reduced to machine round-off. Fig. 32 presents U vectors for the advective equation d 
solution. Frame (b) gives a zoomed in view of the flap nose region. Correctly, the 
vectors are surface normal. Magnitudes are close to unity. 

Fig. 33 gives < 400 ‘error’ histograms. Frame (a) is for the upwind metric 
advective equation. Frame (b) is for the Poisson. The average Frame (a) ‘error’ is 
2.09%. The grid is highly stretched. Therefore, with centered metrics this error rises 
to 10.36%. The Poisson equation has a 1% error. This low Poisson error, relative to 
the advective equation seems to contradict the histogram evidence. This suggests the 
advective d equation has some large ‘errors’ at just a very few points. Clearly, the 
Frame (a) advective equation ‘error’ distribution has the least spread. For the Poisson 
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method (see Frame (b)), it can be seen there is a slight tendency to overpredict d. As 
can be inferred from Table 6, Poisson method lift and drag coefficients are within 
0.05% of those for the NSS procedure. For the advective d equation, the deviation is a 
little greater but not that significant. 



Cl 

Cd 

Advective 

2.815 

0.0582 

NSS 

2.810 

0.0584 

Poisson 

2.811 

0.0584 


Table 6. Case (e) Cl and 
Cl)’s for different d fields. 


Double delta wing (Case (f)) 

For this double delta case the advective d, NNS, and Poisson methods are used. The 
angle of attack is 6^, M = 0.96 and Re = 2.2 x 10^. A key motivation for this case was 
aeroelasticity studies made as part of Boeing’s sonic cruiser development. Fig. 34 
shows initial and severely deformed aeroelasticity calculation surface grids. These use 
1.2 million cells. The meshes are far more severely deformed than for a realistic 
engineering calculation. Extreme plunging generates the deformations. They are 
intended to produce severely deformed cells to strongly test the robustness of 
differential equation based d approaches. Frames (a, b) of Fig. 34 are 3D views. 
Frames (c, d) are 2D y-z plane views. For moving mesh performance studies two 
approximately equi-spaced deflection increments between the Fig. 34 extremes are 
considered. 

Fig. 35 gives advective d equation U component vectors and d contours in three z-y 
planes. Correctly, the Frame (a) vectors are surface normal. The Frame (b) contours 
wrap around the object surface. This suggests their potential grid generation use. Fig. 
36 shows the advective d and Poisson equation residual drops against iteration for a 
stationary mesh solution using the Fig. 34a,c grids. Mesh sequencing with three grid 
levels is used. The abrupt gradient changes generally signify where a finer mesh has 
been moved to. Considering the crude GGS type solver, convergence is adequate, 
requiring around 20 and 40 sec of CPU time for the Poisson and advective d 
equations, respectively. The NSS needs 57 sec. In an aero-elasticity calculation, this 
constitutes around 20% of the time step cost. When solving the Poisson and advective 
d equations over the two equi-spaced deflection increments, just a fraction of a second 
of CPU time is required. This is mostly because of the effectiveness of the spring 
analogy surface grid redistribution procedure (near surface cell relative positions 
hardly change) and also because in a non- stationary grid mode the differential 
equations have a good initial guess for the iterative solution. 

Fig. 37 gives < 400 ‘error’ histograms. Frame (a) is for the upwind metric 
advective equation. Frame (b) is for the Poisson. The average Frame (a) ‘error’ is 
3.27%. With centered metrics this error rises to around 7%. The average Poisson error 
is 2.67%. Encouragingly, the Case (f) lift and drag coefficients for the Poisson method 
are within 0.03 and 0.06% of those for the NSS procedure. 

Table 7 summarises the Poisson and advective equation ‘errors’ for cases (b-f). Prom 
this it is clear that the theoretically exact advective d equation is not always more 
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accurate than the approximate Poisson. This is most probably because the latter is 
more numerically benign. The hyperbolic advective d equation does have the 
advantage that it can be solved in the DES and zonal RANS compatible ME mode 
(also, it can be used for other, non-turbulence related, CED purposes). However, a 
consequence of its hyperbolic form is that upwind metrics are required and also d 
errors are additive (addition occurs in the front propagation direction). The elliptic 
Poisson does not suffer with this additive d error trait. As discussed in Section 2, the 
advective equation can be cast in a Hamilton-Jacobi form. This gives the potential to 
control convex and concave feature d distributions, and so remedy turbulence model 
defects. Although the Poisson method offers no explicit convex/concave feature d 
control, its natural traits in such regions appear helpful. It is difficult to conclude 
which method is best. This is especially the case since their use is purpose dependent. 
For example, the Poisson approach is totally unsuitable for evaluating far-field 
computational interfaces on overset meshes as needed for the approach of Nakahashi 
and Togashi (2000). It is also worth remembering that ‘error’ values given here are 
relative to the NSS approach, and that this cannot be regarded as providing a perfect d 
field. 


Case 

% ‘Error’ 


Advective 

Poisson 

(b) 

4.8 

17 

(c) 

0.5 

3 

(d) 

3.1 

3.1 

(e) 

2.1 

1 

(f) 

3.3 

2.7 


Table 7. Case (b-f) d error summary. 


17 




5. CONCLUSIONS 

Advection and advection-diffusion forms of the Eikonal and Hamilton- Jacobi wall 
distance {d) equations have been presented. These have been studied because of their 
relatively easy potential implementation in industrial CFD solvers. Mostly, use of the 
NASA CFL3D CFD program to solve the advective -based d forms was explored. The 
advection based d equations were found to have robust convergence. Geometries 
studied included single and two-element airfoils, wing body and double delta 
configurations, along with a complex electronics system. It was shown that for 
Fikonal accuracy, offset metric differences are required. A simple, approximate, 
Poisson-equation-based distance approach was found effective and, since it did not 
require offset metric evaluations, easiest to implement. The sensitivity of flow 
solutions to d assumptions was explored. Generally, results were not greatly affected 
by wall distance traits. However, the use of, for example, nearest surface search 
distances was found to reduce downstream trailing edge intensities by around 4.5%. 
Results suggest that for stationary meshes the advective d equation has about twice 
the computational cost of the Fikonal, but for moving meshes this difference greatly 
decreases. However, CFLSD’s hybrid TFFspring analogy grid redistribution scheme 
maintains near wall relative grid positions well. Consequently, for non-stationary 
grids, differential equation based d approaches were not greatly challenged. 
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Fig. 2. Case (a) grids: (a) mildly distorted and (b) strongly distorted. 
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Fig. 3. Predicted against actual distances (y) at x=0.5 for flat plate: 
(a) strongly stretched grid and (b) strongly distorted 
( exact solution, O offset metrics, A centered metrics). 
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Fig. 4. Velocity data for strongly distorted grid: (a) lUI and (b) U vectors. 
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Fig. 7. Case (b) (Electronic system) mid-plane grid structures: 
(a) x-y plane and (b) z-y plane. 
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Fig. 11. Case (b) contours of d in the x-y plane for the following approaches: (a) 
Search (NNS); (b) Poisson; (c) Eikonal; (d) advective; and (e) advection-diffusion. 
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Fig. 12. Case (a), influence of Laplacian: (a) Region A and (b) Region B. 
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Fig. 13. Case (b) (electronics system) RANS/LES interface for a DES solution using 
the advective d equation with the ME approach. 
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Fig. 17. Case (c) (NACA 4412) d contours: (a) NSS; (b) NNS; (c) advective equation; 
(d) Poisson and (e) overset grid Poisson. 
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Fig 18. Case (c) (NACA 4412) < 400 Error histogram: (a) advective d equation and 

(b) Poisson. 
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Fig. 19. Case (c) (NACA 4412) Poisson d ‘error’ contours. 
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Fig. 26. Case (c) (NACA0012) pitching wing results: (a) lift coefficient variations 

and (b) moment coefficient variations. (O measurements, non-deforming grid, _ _ 

deforming grid with NSS, . . . . deforming grid with Poisson) 







Fig. 27. Case (d) (wing body) grid. 
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Fig. 30. Case (e) (wing flap) comparison of d contours: (a) search (NSS); (b) 

advective; and (c) Poisson. 










: (a) advective equation and (b) Poisson. 






g. 34. Case (f) (double ( 
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Fig. 35. Case (f) (double delta) advective d equation results 
in three z-y planes: (a) U vectors and (b) d contours. 
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Fig. 37. Case (f) (double delta) y < 400 J error histogram: (a) advective equation and 

(b) Poisson. 
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